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5^ I (February 1, 2008) 

O^ \ A detailed study of the Frohlich polaron model is performed on the basis of diagrammatic quan- 

tum Monte Carlo method W. The method is further developed both quantitatively (performance) 
and qualitatively (new estimators), and is enhanced by spectral analysis of the polaron Green func- 
tion, within a novel approach. We present the best up to date results for the binding energy, and for 
the first time make available precise data for the effective mass, including the region of intermediate 

^+ ' and strong couplings. We look at the structure of the polaron cloud and answer such questions as the 

average number of phonons in the cloud and their number/momentum distribution. The spectral 
analysis reveals non-trivial structure of the spectral density at intermediate and large coupling: the 
spectral continuum features pronounced peaks that we attribute to unstable excited states of the 



flj ■ polaron. 

PACS numbers: 71.38.+i, 02.70.Lq, 05.20.-y 
-(— > 



I. INTRODUCTION 



The polaron problem (for an introduction see Ref. 0) has originally emerged in the solid state physics as a problem 

of electron moving in a (dielectric) medium. It became clear, however, that this problem is of essential general-physical 

interest, as a model of a quantum object strongly coupled to an environment. Starting from the work of Landau [pj, 

2 . the polaron problem has been attracting a permanent attention, serving as a testing ground of new non-perturbative 

Q ^ methods. 

The most popular model in the polaron problem is the so-called Frohlich Hamiltonian describing an electron coupled 
to non-dispersive (optical) phonons of a dielectric medium via its polarization (Plank's constant and electron mass 
J> . are set equal to unity) ||] : 

in 

(N : iJ = i/e + ffph + ffe-ph , (1.1) 

O 

0^■ H, = Y,\<^W. (1-2) 

0^■ 



^Ph-E'^^^q^' (1.3) 

q 



S ■ ^-ph = Y. ^(q) (^q - ^-q) «k-q«k , (1-4) 

k,q 



^ ■ / ^ \l/2 1 

'ij ■ Viq) = i (2\/2Q7r) - . (1.5) 

In Eqs. ( |l.l| - [l.5| ), a^ and h„ are the annihilation operators for the electron with momentum k and for the phonon with 
momentum q, respectively, ojq = wq is the q-independcnt phonon frequency, which can be set equal to unity without 
loss of generality, a is a dimensionless coupling constant. 

Despite a lot of work addressed to the Frohlich Hamiltonian, the model is still far from being completely understood. 
In the most interesting region - at intermediate and large values of a, almost all available treatments are of variational 
character. Hence, these treatments, even if consistent with each other, can not guarantee quantitative and qualitative 
reliability of the results. Moreover, some treatments are known to be in qualitative disagreement with the others. 
As a characteristic example, note that certain approaches suggest that the polaron states at small and large a's are 



of qualitatively different nature, and there should occur a sort of phase transition in the parameter a [|5|-|ll| (see, 
however, discussion in Ref. |1^). The study of such important issues as electron ^-factor and the structure of the 
polaronic cloud was also restricted to the perturbation theory and variational treatments at small momenta |13|-p^ . It 
remained unclear what are the limits of applicability of these results, and whether they correctly describe the physics 
of polarons in the most important range of intermediate a. 

Recently, a method of diagrammatic quantum Monte Carlo (MC) was developed, which is very efficient for the 
polaron-like problems B. The method allows direct simulation of entities specified in terms of (positive definite) 
diagrammatic expansions. In Ref. ||I| the polaron Green function was simulated, and the results were used to extract 
the polaron spectrum. 

In the present paper, we employ the diagrammatic Monte Carlo scheme of Ref. [0 for a detailed study of the Frohlich 
model. We significantly enhanced the original scheme by (i) introducing 7V-phonon Green functions (with 2N external 
phonon lines), which are simulated in one and the same MC process with the ordinary (0-phonon) Green function, (ii) 
developing a powerful procedure of spectral analysis of the Green function. The A^-phonon Green functions allow us to 
consider the structure of the phonon cloud and facilitate obtaining polaron parameters at large a, where the polaron 
is essentially a many-phonon object. In particular, direct estimators for the energy, effective mass, group velocity, 
and Z-factors can be constructed. The spectral analysis of the Green function gives the most complete information 
about the polaron, including the possibility to reveal stable and metastable excited states, if any. 

The paper is organized as follows. In Section || we introduce the set of Green functions, describe the corresponding 



diagrammatic series, and discuss how they are related to the polaron parameters. In Section [II we describe qual- 
itatively the Monte Carlo procedure [the quantitative discussion of the updates is given in the Appendix A]. The 
calculated properties of the polaron (energy, effective mass, structure of the polaronic cloud, ets.) are presented in 
Section IV. In Section ^ we analyze excited states of the polaron by restoring the spectral density of the Green 
function by a novel method of the spectral analysis [detailed description of the method is presented in Appendix B] . 

II. GREEN FUNCTIONS AND DIAGRAMS 

In this section we introduce basic entities and establish their relations, which will be utilized in the rest of the 
paper. 

We start with the standard Green function of the polaron in the momentum (k) - imaginary-time (r) representation: 

G(k, r) = (vac| a^{T)ai{Q) |vac) , r > , (2.1) 

«kM = e"^a^e-"^ . (2.2) 

Here |vac) is the vacuum state. 

The physical information that G(k, r) contains is clear from the expansion 

G(k,r) = Y. l(H«|.|vac)|2e-(^^('')-^'')- , (2.3) 

where {\y)} is a complete set of eigenstates of the Hamiltonian H in the sector of given k, i.e. H |i^(k)) = Ey{k) |i^(k)), 



H |vac) — E{) |vac). Since in our model Eq — 0, we omit it below. Rewriting Eq. (|2.3| ) as 



oo 



G(k,r) = / dL^5k(c^)e-'"^ (2.4) 







g^iuo) = Y,5{uj- E,{\^)) |(t^|aj,|vac)|2 , (2.5) 

one defines the spectral function 5k (^) which has poles (sharp peaks) at frequencies corresponding to stable 
(metastable) particle-like states. Hence, if at a given k there exists a stable polaron with the energy i?(k), the 
spectral function reads 

gk(a;) = Z^^) b(i^ - i?(k)) -I- . . . , (2.6) 

where 



/ 

/ 

/ 


/ 

/ 
1 


/ 

/ 


' \ 
\ 

\ 


\ 


\ 
\ 


\ 

\ 
1 


/ 


- 


- 


- 


\ 


; 4 


» 


— — • — 


— * 




-♦— — 


*— 






__ 


__ 


— * — c 



FIG. 1. A 0-phonon diagram. 
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FIG. 2. A 2-phonon diagram. Note, that disconnected phonon lines always appear in pairs (labels A and B), both lines in a 
pair having the same momentum. 



Kk) 



polaron (k) | free electron (k) ) p . (2.7) 



Moreover, if the polaron state is the ground state, its energy and Z-factor are "projected out" by the Green function 
behavior at long times: 

G(k,r>Wo"') ^ Z^^^ e-^'^^^^ . (2.8) 



Along with the standard polaron Green function (2.1), it is reasonable to introduce the A^-phonon Green function 

N 

G^(k,r;qi,...,qA,) = (vac| &q„(r) • • • 6^^ (r) ap(r)at,(0) 6^^(0) • • • &^„(0)|vac) , p = k - ^ q, . (2.9) 



Relations (2.3-2.8) are readily generalized to the case of A^-phonon Green function. In particular, the A^-phonon 



Z-factor for the stable (groundstate) polaron with momentum k, 

Zn {^1-, ■ • • I qjv) = I ( polaron (k) | free electron (p) + free phonons (qi, . . . , qjv) ) p (2-10) 



[the momentum p is defined as in Eq. (2.9)], is given by 



GA,(k,r»Cc;o-i;qi,...,qjv) ^ ^(.^^(qi, . . . , q^) e^^^*')" . (2.11) 

Our MC procedure of simulating Green functions will utilize a standard diagrammatic expansion - Matsubara 
technique at T = 0. The diagrams (see Figs, ll 0) are built of the following elements: (i) free-electron propagator 
(solid line) 



G("'(p,r2-ri) = exph^(T2-ri)] , (2.12) 

(ii) phonon propagator (dashed line) 

D{q, T2 - n) == exp[-ujq{T2 - Tl)] , (2.13) 

(iii) vertex factor V{q) ascribed to the vertex formed by a phonon propagator (with the momentum q) and two 



adjacent electron propagators. External lines of diagrams arise from the operators standing in Eqs. (2.1), ( p.9[ ) and 
their times and momenta are defined accordingly. Momenta of the internal lines are free up to the momentum 
conservation constraint at each vertex. To obey this constraint, we choose phonon momenta to be free, fixing thus 
the momenta of the electron lines. Times ascribed to the line ends are subject to the chronologization constraint: 
The time of the left end is 0, the time of the right diagram end is r, times of the vertexes must increase along the 
global electron line, directed from to r. Integration over all free parameters of the diagrams is assumed, within the 
domains consistent with the constraints. 

Upon the formulation of the diagrammatic rules, it is reasonable to introduce the irreducible A'^-phonon Green 
function, Gn, which consists only of irreducible diagrams. [By irreducible diagrams we understand those ones that 
do not contain phonon lines decoupled from the electron.] ^From this definition it is clear that the reducible part of 
the iV-phonon Green function, Gn — Gn, is a sum of products of irreducible functions Gn' , N' < N, and free phonon 
propagators. Therefore, the reducible part does not contribute to the large-r asymptoti cs of Gn [because of extra 
factors e~'^°^ coming from the disconnected phonon propagators], and, in particular, Eq. ( p.ll ) holds true for Gn as 
well. 

It is useful! then to consider the following function: 

oo ^ 

P(k, r) = G(k, t) + Y, dqi-- dqN Gjv(k, r; qi, . . . , qjv) (2.14) 

N=l'' 

[Note, that if Gn -^ Gn , this expression would be singular at r ^ because of the divergence of the integrals for 
disconnected phonon propagators]. The function P is readily calculated by our MC procedure (see next sections) and. 



according to (2.8), (2.11), satisfies the relation 

P(k,r»co-i) ^ e-^('')^. (2.15) 

Here we took into account the completeness relation 

oo „ 
^o"^ + E dm---dqN Z^j^\q,, . . . ,qAr) = 1 . (2.16) 

N=l-^ 



Eqs. ( 2.15| ), ( ^.16 ) imply that the lowest level -E(k) is non-degenerate. Otherwise one should introduce degeneracy 



factors to the right-hand sides of the equations. 

III. DIAGRAMMATIC QUANTUM MONTE CARLO 

In the previous work [yjjlj] it has been shown how to sum convergent (and arbitrary otherwise) diagrammatic series 
numerically without systematic errors. In this section we outline the basic numerical procedure of evaluating series 
for various Green functions described in the previous section [the details being discussed in the Appendix A], and 
introduce a number of estimators, that render the evaluation process significantly more efficient. 

Suppose that we are interested in the function Q{{y}), which depends on a set of variables {y}, and which is given 
in terms of a series of integrals with an ever increasing number of integration variables: 

oo „ 

Qi{y}) = ^ ^ / da;i---da;ml?m(^m,{2/},a:i,...,a;™) . (3.1) 

Here ^m indexes different terms/diagrams of the same order m. The term ttt, = is understood as a certain function 
of {y}. Both external {y} and internal {xi} variables are allowed to be either continuous or discrete; in the latter case 
integrals are understood as sums. Diagrammatic MC process is a numeric procedure based on the Metropolis principle 



| jl8[ that samples various diagrams in the parameter space ({y}, m,^m, {x}m) and collects statistics for Q{{y}) in such 
a way that the final result converges to the exact answer. The process has very much in common with the Monte 
Carlo simulation of a distribution given by a multi-dimensional integral. Nevertheless, there is an essential difference 
associated with the fact that integration multiplicity in the expansion Eq. (pj) is varying. 

Summing the series for Q{{y}) is the process of sequential stochastic generation of diagrams described by functions 
Vm, Eq. (pj). The MC process consists of a number of elementary updates falling into two qualitatively different 
classes: (I) those which do not change the type of the diagram (change the values of arguments in "Dm, but not 
the function itself), and (II) those which do change the diagram order. The set of elementary updates and their 
implementation is problem-specific; the only necessary requirements are ergodicity, i.e., given two arbitrary diagrams 
it takes finite number of updates to transform one to another, and detailed balance, i.e., diagrams contribute to the 
statistics according to ratio of their P-functions. In the Appendix A we describe a set of updates which we find to 
work very efficiently. However, considering enormous freedom in constructing various updates, we have little doubts 
that they may not be further improved. 

Though Green function contains complete information associated with the polaron spectrum, and the accumulation 
of the histogram for G is straightforward, it is reasonable to introduce certain direct estimators, including one for 
the Green function itself. These estimators substantially enhance the accuracy of calculations and/or allow collecting 
more information during one MC run (by "spreading" the data to the different values of the external parameters). 

A. Estimators for effective mass, group velocity, and energy 

We start with the family of estimators that are constructed in accordance with the following standard MC rule. 
Suppose we have some quantity A specified by the diagrammatic expansion 



A = J2 T^i^^ ' (3-2) 



where T)], s are the diagrams for A which are parametrized by internal variables denoted by the unified index i/, 
and the summation over v is understood as the summation over discrete variables and integration over the continuous 
ones. Suppose next that all I?i, s are positive definite and we have a MC process of generation of random v with 
probability density given by I?}/ . Then, if some quantity B is specified by a (similar) diagrammatic expansion 

B = Y^ Vi^) , (3.3) 

the estimator for the ratio B /A is given by 



Here mca{v} means the set of i^'s generated during the MC run. Commonly, the quantities A and B in Eq. (3.4) 
are, respectively, the partition function and an observable. We, however, will use this relation in a somewhat different 
context. For one thing, in our case A and B can correspond to one and the same Green function, but at different 
values of external parameters (say, momentum or coupling constant). This way we are able to obtain results for a 
number of different values of the external parameters from a MC process for just one fixed set of parameters. We can 
also directly calculate derivatives with respect to the external parameters. To this end we should analytically take 



the corresponding limit from both sides of Eq. (3.4) 



Let us obtain estimators for the effective mass, ?n*, and group velocity, v(k) := 9-B(k)/9k. First we note that 
(r ^ oo, A ^ 0) 



P(k + Ae,r) / exp(-A2T/2m*) , fc = 0, 



P(k, r) \ exp(-Aev(k)T) , k^Q, 



(3.6) 



where e is a unit vector. Considering the denominator and the numerator of the left-hand side of Eq. (3.6) as A and 
B (respectively), we can take advantage of Eqs. (3.4) and (3.5). The function Q is given by 



Q = ncxp{-kp.- + A6)2 - P,'](Ar),} . 



(3.7) 



Here j numerates free-electron propagators of a given diagram for P(k, t); Pj is the momentum corresponding to 
the propagator j, and (Ar)j is the length of this propagator. Eq. (^J) immediately follows from the fact that the 
series for P(k -f- Ae, r) can be obtained from the series for P(k, r) by adding the momentum Ae to all free-electron 
propagators. As we are interested only in the limit A — > 0, we can expand Eq. (3.7) in powers of A: 



Q - 1 - Ar(ep) - y r + ^r^iepf + 0{X') 



where p is the mean electronic momentum of the given diagram 



(3.8) 



(3.9) 



Comparing Eq. (3.8) with the corresponding expansions of the right-hand side of Eq. ( |3.6| ), we arrive to the following 
estimators 



'MC 



v(k) (t -^ oo) 



(3.10) 



1-3 ((pH 



MC 



1 



(t ^ oo) , 



(3.11) 



where ( . . . )mc nieans MC averaging in accordance with Eq. (3.4). 

A special care should be taken for treating the time of the Green funct ion as a n external parameter of the diagrams 
(in the sense adopted in this section). The problem is that relations (3.4-3^51) imply that the internal parameters 
of diagrams VI and Vl have one and the same domain of definition, otherwise the ratio (3.5) is not correctly 



defined. Meanwhile, the domain of internal times of diagrams directly depends on the external time. To circumvent 
this problem, one can introduce scaled internal times by simple relation r^ = rfi , where r is the external time (length 
of diagram in time), t^ is an internal time variable (position in time of an electron-phonon vertex), and f^ g [0, 1] is 
the corresponding scaled time variable with the domain of definition independent of t. 

Now it is easy to obtain a direct estimator for the polaron energy. To this end we start from the relation 



PjK (1 + A)t) 
P(k,T) 



,-AB(k)r 



(r -^ oo) 



(3.12) 



and proceed analogously to Eqs. (3.6-[3.11). In this case for the function Q we have 



Q = (1 + A)^ m^^pt- 



-A|(Ar),] 



J|cxp[-AtJo(AT) 



(3.13) 



where indexes j and s stand for the electron and phonon propagators, respectively, and N is the number of integrations 
over times (or, equivalently, number of interaction vertexes) in a given diagram. Then, in the limit A — > 0, we expand 
the right-hand sides of Eqs. (3.12- 3^1^ ) up to terms proportional to A, and in accordance with Eqs. (3.4-3.5) arrive to 
the estimator 




(Mj + J2 '^o^^^) 



N 



E(k) (r ^ oo) 



(3.14) 



MC 



B. Reweighting 

We will also employ the reweighting technique ||l^, which allows one to utilize the statistics being generated for 
some given set of external variables, £,, for calculations at a different set, ^'. In terms of the diagrammatic Monte 
Carlo, this technique is based on the relation 

E Q'^io - E iil|Q-(^')' (3.15) 

where Q^, is any quantity summed over MC statistics. [We omitted superscript A at V^, since it is not relevant here.] 
The relation (3.15|) follows from the fact that the MC statistics for the set S,' involves the same (in the sense of 



structure and the values of internal parameters) diagrams as the statistics for the set ^. The difference is only due to 
a different probability to generate a diagram with the set ^ rather than ^'. This difference can be taken into account 



analytically by the corresponding ratio, which immediately leads to (3.15). 

In our case, typical external parameters are the interaction constant and the polaron momentum: ^ — (a, k). The 
corresponding ratio of the diagrams is 

^^ - (^) n-P{-[(k' - k)^ + 2p,(k' - k)](Ar),/2} = (^j ,-[(k'-k)V.+P(k'-.)]. . (3.16) 

This relation allows to get many points at different a' and k', at no extra cost in CPU time, while performing MC at 
a given set of a and k (cf. Ref. [|l9||). 

C. Exact estimator for Green's function 

Calculation of the Green function by means of histogram, though is simple and natural, involves an apparent 
shortcoming, associated with the finite width of the histogram cell. There is always a competition between the 
decreasing systematic error by making the size of the cell smaller, and the increasing statistical accuracy, which 
requires increasing the size of the histogram cell. This problem can be solved by introducing an exact (free of 
systematic errors) estimator for the Green function, as follows from a generic consideration presented below. 

Given some function A(^o) of an external variable/set of variables ^Oj specified with the (positive definite) diagram- 
matic expansion 

A{^o) = Y.V4^o)^ /d^E^'^(^)'^(^-^")' (3.17) 

[considering a general case, we do not assume that the domain of definition of i/ is independent of ^] and having 
arranged a MC process of generating configurations {i^, C} with the probability density proportional to 'D^{^), we 
would like to construct an estimator a^g the average of which over the MC process gives (up to a global normalization 
factor) the function A{^o)- Let us look for a^g in the following form 

a,.M = I f )^^(^°)/^^(^)' \^ ^ r„ and P.(0 ^ , 
5o\ '>/ 1^ Q ^ otherwise. 

Here Fq is some finite domain in the space of variable ^ including the point ^o, <l{i^) is some function to be defined 
later. [We adopt a convenient and consistent with the MC procedure convention that 'Di,{(^) = 0, if ^ is out of the 
range of definition of the corresponding diagram.] /^From Eq. (pTq) we have 

(a«o>MC ^ ^E I diaUv,OV,{i) = C^ g(^)I?.(^o) / d^ (3.19) 



where 



C' = Y. J d^V.iO (3.20) 



is the normalization factor for the distribution of the random pairs {v, £,) induced by the series ( [3.17 ). ^'^From Eq. (3.1£) 
it is seen that if we choose 

q-'H - / d^ (3.21) 



note that, according to ( 3.18 ), the definition of q{i') is relevant only when I'iy(Co) 7^ 0, and that q "^{v) ^ 0, since at 



least small neighborhood of the point ^0 contributes to the integral] , then 

(«eo)MC ^ ^^(^0). (3.22) 

The particular form of the estimator for the Green function is readily obtained by identifying ^ with r, and noting 
that the ratio of diagrams standing in Eq. (3.1^) is given in this case by Eq. ( 3.13| ). As the domain of definition of 



any diagram with respect to r is independent of the diagram structure: r e [0, 00], the factor q in Eq. (3.18) is simply 
proportional to the inverse size of the interval Fq. The choice of Fq for each particular tq is arbitrary, being a matter 
of taste and convenience. 

D. Improved estimators for phonon statistics 

Collecting statistics for the phonon cloud can be significantly improved by a trick described below. 

We start with noting that in the case of r -^ 00, which is relevant to the ground-state properties, the set of all 
iV-phonon diagrams possesses a certain symmetry. To reveal this symmetry, we transform diagrams to the circular 
representation by the rules illustrated in Fig. y. The transformation involves "gluing" the outer ends of the electron 
line, as well as the outer ends of pairs of (corresponding to each othe r) ex t ernal phonon lines. It is easy to check that 



the procedure is consistent with the definitions of propagators, Eqs. ( 2.12 , 2.13 ). In the limit of large r, which we are 
interested in, the probability to find a phonon propagator with length > r/2 is vanishingly small. That is why we 
may omit orientation, understanding the time length of phonon propagator as the length of the smallest arc between 
its ends. With the same accuracy, in the limit r — > 00, we may ignore the pairs of phonon lines like the one shown in 
Fig. 0(c) , which do not have unambiguous circular counterparts. 

It is worth noting that circular diagrams naturally occur in a finite-temperature technique (cf., e.g., Ref. pC|), 
where the circumference r has a meaning of inverse temperature. 

Within the circular representation the symmetry necessary for constructing improved estimators is clearly seen. 
Indeed, a circular diagram represents a whole class of plain diagrams, due to its independence of the position of the 
point corresponding to the ends of the plain diagram. 

Thus, having generated a certain plain diagram and associating it with the corresponding circular diagram, one 
effectively produces a whole class of diagrams to be included into the statistics. 

In practice, the procedure is as follows. Let index i = 1, . . . , N^ label electron propagators in the circular diagram; 
we thus split the diagram into N,, pieces each having duration Ar^ in time, X)i=i ^i ~ '^- When the circular diagram 
is cut anywhere on the interval Ar^ we obtain a contribution to GjVi(k, riq^ , . . . iQtv )' where iV^ is the number of 
phonon propagators which are cut along with the electron propagator on this interval, and {qi} are their momenta. 
An estimator for the integrated N-phonon Z-factor is found then to be 

Z^^^^ ... II dq,Z^^\qi, . . . , q^) = / ^ — <5^.,Ar ) , (3.23) 

i.e., due to time invariance of the circular representation each interval contributes to the statistics according to its 
duration in time. The one-phonon distribution function within the N-phonon states manifold is given by 

/» w -. TV » » TV 

• • • / n d(\]Z'~-^\(\, q2, ■ • ■ , qjv) = ^ ^ / . . . / n dq.jZ'^^\(ii, . . . , q; = q, . . . , qw 




MC 



Obviously, Z)^ = J dqFj^ (q) . Summing over all N we obtain the one-phonon distribution function 



(a) 



(b) 



(c) 




(a') 




(b') 



FIG. 3. Correspondence between the standard and circular representations of diagrams. 
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FIG. 4. Bottom of the polaron band Eo as a function of a. The error bars are much smaller than the point size. 



i^(k)(q) 



N=l 



(3.25) 



IV. NUMERIC RESULTS 



A. Groundstate energy and effective mass 



In Fig. we present our results for the bottom of the band Eq as a function of a in a wide region of coupling 
strengths. The data are compared with Feynman's variational treatment p3| to demonstrate the remarkable accuracy 
of the Feynman's approach to the polaron energy. We thus conclude that just a simple extrapolation between 
the second-order perturbative result E^ — —a — 1.26(q;/10)^ and Feynman's strong-coupling variational estimate 
Eq = —o? I'i'K — 2.83 [such an extrapolation is very close to Feynman's variational treatment in the whole range of 
a's] yields quite satisfactory approximation for Eq{(x). 

In Fig. |5| accurate data for the effective mass are presented up to m, ~ 1000 [For larger m,'s the Frohlich model is 
not realistic anyway since the phonon dispersion becomes relevant] . At a < 9 the statistics were collected at integer 
values of a with reweighting (in accordance with the procedure described in the previous section) to corresponding 
finite intervals (see the plot). At a > 9 the reweighting procedure proved to be ineffective. 

Let us compare the data for m, with the small- and strong-coupling analytic results. At small a's, the formula 
TO* = (1 — a/6)~^, known to coincide with the perturbation expansion up to the second order |}4|, |l6yi4| , works well 
up to a « 2. In contrast to the case of i?o, the strong-coupling limit [|l| to.* = a'*/48 drastically overestimates the 
effective mass in the whole range of physically interesting a's. Feynman's variational technique works better, but still 
with a considerable deviation (up to 50% ) in the region 5 < a < 10. 

Almost the same degree of accuracy gives variational treatment by Feranchuk, Fisher, and Komarov |ll[ . An 
important point about this treatment, however, is that it suggests a phase transition from "light" to "heavy" polaron 
at a close to 7.5, which should lead to an especially rapid increase of ?7i* just after a = 7.5. In this connection, we note 
that our curve to* (a) is essentially smooth and does not suggest any sort of phase transition or sharp cross-over. [The 
results of more deep numeric study of the possibility of the phase transition are presented in the next (sub)sections.] 
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FIG. 5. Effective mass as a function of coupling parameter. Our MC data (circles interpolated by solid line) are compared 
with perturbation theory and strong-coupling-limit results (dashed lines), Feynman's approach (squares), and Foranchuk et. 
al. variational approach (dimonds). 
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FIG. 6. The bare-electron Zq -factor for the ground state as a function of the couphng strengt h; fi lled circles are the MC 
data (with error bars smaller than point size), and the solid line is the perturbation theory result ([4.l|). 



B. Structure of polaronic cloud 



(k) 

In this subsection we present our results for the electron Z-factor, or Zg , for k = as a function of the coupling 
strength in the region of small and intermediate a (for a > 10, the bare electron weight in the polaron ground state 
becomes vanishingly small), and for a = 1 as a function of momentum up to the end-point. We study also, how the 
distribution of phonons in the ground state, Zj^ , and the average number of phonons, N = X)w=i ^^n ' evolve with 
a. Finally, we show how the physics of the end-point is seen in the transformation of the one-particle distribution 
function -Fklq), Eq. ( t?.2[ ), and how it allows to identify the relevant self-energy diagram. 

For small a and k = the leading behavior is readily obtained from the perturbation theory 



Fr(q)=^" 



smt 



47r2 [g2/2 + l]2 



dq d9 dip 



Ao) 



a/2 



Zr = l-a/2 



(4.1) 



We have verified that for a < 1 the perturbative results (4.1) are describing the data rather accurately [see also 
Figs. 0, 0, and dashed curves (connecting filled circles) in Fig. |11|] . 

The data in Fig. M make it clear that perturbation theory may not be trusted for a > 1 when the bare-electron 
state in the polaron wave function is no longer the dominant contribution, e.g., Zq {a = 3) < 0.2. The bare electron 
Z-factor vanishes rather rapidly for a > 3 (the dependence ZqIo) is faster than exponential) and becomes < 10~^ 
for a > 10. We do not attempt to fit the data to the particular functional dependence since we believe that in the 
interval 3 < a < 10 the polaron state undergoes a smooth transformation between weak and strong coupling limits. 

In Fig. M we show the distribution of multi-phonon states in the polaron cloud at A; = 0. We see how it gradually 
evolves from the perturbation theory case into the strong-coupling regime. For a > 10 the data may be fit to the 
Gaussian distribution, but at smaller values of a the distribution is essentially asymmetric - it decays faster for N > N 
than for N < N. We note, that even for alpha = 17 the phonon number distribution is rather large which means 
that the polaronic cloud is essentially a superposition of states with different N. The effects resulting from this fact 
are outside the scope of the variational ^^-theory and may, for example, account for the considerable deviation of the 
effective mass discussed above. 
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FIG. 7. Partial contributions of A'^-phonon states to the polaron ground state for various values of a. Error bars are shown, 
but are typically smaller than the point size. (The dashed lines are to guide an eye.) 
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FIG. 8. The average number of phonons in the polaron ground state as a function of a. Filled circles are the MC data (error 
bars are smaller than the point size), the dashed line is the perturbation theory result (i.l), and the solid line is the parabolic 
fit for the strong coupling limit. 



It is well known that polaron models often support the self-trapping phenomenon, when the ground state changes 
in a relatively narrow interval of parameters from light- to heavy-mass state with a sharp increase in the number 
of phonons contributing to the polaronic cloud. The same phenomenon was advocated for the Frohlich model by 
a number of authors P-ni|, including the statement that more than one stable polaron state exists in the region 
of intermediate a. Clearly, if there were a sharp transformation of the polaronic ground state, it would have been 
immediately seen in the phonon statistics. Such a transformation might be not visible on the energy plot if the 
hybridization matrix element between the two competing states is not small, and their energy derivatives, dE/da, 
are close to each other at the point of crossing. However, if we are to speak about different polaronic states, then 
(almost by definition) their structure has to undergo an abrupt change with a. Our data on the bare electron Z-factor 
Zq and phonon distribution functions Z^ are evolving smoothly with a and thus prove continuous formation of the 
self-trapped state. 

To further support this conclusion, we plot in Fig. Hthe dependence N vs a. The crossover between the perturbative 
result N ~ a/2 of Ref. p6| and the strong-coupling limit, where N ^ 0.22a^ , demonstrates no sign of the level crossing 
picture. As a side remark we note that the result of Ref. p6| which predicted that perturbation theory for N works 
well in the intermediate range 1 < a < 6 is not true. In fact, this law breaks down along with the perturbation theory. 

Consider now the evolution of the polaronic cloud with momentum as we approach the end point. [The dispersion 



curve E(k) featuring the end point at momentum kc of the form E(k — > kc) = Eq - 



-UJQ' 



{k — kcY /2mc was calculated 



in Ref. ||l|.] Although the polaronic state is stable for E{k) ~ E(0) < luq the bare electron weight vanishes as fc ^ k^ 
see Fig. H. From this plot we estimate kc{a = 1) w 1.83. This figure also makes it clear to what degree the earlier 
result that Zq is momentum-independent |22(| works. 

With the numerical tools at hand it is possible to "visualize" the physics of the end point. In accordance with 
the generic Pitaevskii theory of the end point p3[ , in the vicinity of the end point the polaron can be considered 
as a (weakly) bound state of phonon, carrying almost all the momentum of the state, and a polaron with almost 
zero momentum. This physics is transparent from the comparison between the statistics of A'^-phonon states in the 
ground state and at A: — *■ fcc shown in Fig. nu. Evidently, the two curves can be matched by shifting the ground- 



state distribution by one, i.e. Z 



(fee) 

N 






7(fe=1.79) 



(the maximum of the Zj^~'"'^' curve is a little depleted because of the 
remaining finite weight Zq since k < kc). It means that the polaronic state near the end point is a superposition of 
bound states of the phonon with momentum around kc and a polaron at the band bottom. 

Since Fig. |lO| does not tell us explicitly what are the parameters of the extra phonon present in the polaronic cloud 
at A: -^ kc, we plot in Figs. 11 12| normalized distribution functions of phonon momenta (in |g| and in the angle 
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FIG. 9. The bare-electron Z-factor as a function of the polaron momentum up to the end point, 
between qk). It is obvious from these figures that the extra phonon momentum is concentrated around k. 

V. SPECTRAL ANALYSIS 

The spectral function g-k{oj) ( ^.5D provides important information about the system since it has poles (sharp peaks) 
at frequencies corresponding to stable (metastable) particle-like states. Besides, since the probability of absorption of 
a free electron with the momentum k into a polaron state is proportional to the spectral function, the latter can be 
measured experimentally by angle-resolved inverse photoemission spectroscopy. 

In spite of many elaborated treatments of the properties of the polaron, the knowledge about high-energy part of 
the polaron spectrum is mostly limited by attempts to calculate the spectral density either by perturbation theory 
approaches or at strong coupling limit p4| . As both the Green function asymptotic behavior and the machinery of 
estimators provides information about ground state properties only, the spectral density is indispensable for the study 
of excited states of the system. 



The spectral analysis, i.e. solving the equation (2.4), was performed by a novel method [detailed description of the 
method and testing examples are presented in Appendix II] . The most important features of the method are that it 
avoids distortion of equation by nonlinear terms and does not suffer from systematic errors caused by preassigned 
discretisation of the w-space. 

To perform a joint check of the diagrammatic Monte Carlo approach and the method of spectral analysis, we 
compared the spectral densities obtained by our numeric calculations and by perturbation theory for zero temperature. 
The analytic expression for the high-energy part (a; > 0) of the spectral density could be obtained for the arbitrary 
interaction potential V(\ q |) which depends on the modulus | q | of the phonon momentum. For zero polaron 
momentum, k = 0, the the imaginary part of the linear in a self-energy part S(0, w > 0) is 



lml](0,cj) 



\/27r 



V^T^T I v{y/2{uj-i)) p e{uj - 1). 



(5.1) 



(Here 9 is the theta-function. 
to a terms one gets 



Then, using the relation go{uj) — — ImGk=o('-^)/7r and keeping only linear with respect 



goiuj > 0) = 



1 yiJ^^ 

V27r2 LU^ 



\V{VW^^)fO{cu-l). 



(5.2) 



The expression for the low-energy part {lo < 0) of the spectral density depends on the specific form of the interaction 
potential and we consider the perturbation theory result for the short-range interaction 
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FIG. 10. Partial weights Zn of the A'^-phonon states in the structure of the polaronic state for a = 1 at the band bottom 
and near the end point. (The dashed hues are to guide an eye. 

F(q)r 




FIG. 11. Phonon distribution functions in g-modulus for the ground state (filled circles) and close to the end-point for 
k — 1.79 (open circles). The momentum fcc is indicated by a bar at the q-axis. (The lines are just linear interpolations between 
accurate numeric points.) 
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FIG. 12. Phonon distribution functions in the angle between the vectors q and k close to the end-point for k = 1.79. 
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FIG. 13. The comparison of the numeric results (solid lines) and the perturbation theory curves (dashed lines) for the spectral 
density of Frohhch model with a = 0.05 (upper panel) and the short-range interaction model with a — 0.05 and k = 1 (lower 
panel) . 



V{\ q I) = i (2V20 



1/2 



1 



vV+^ 



which reduces to the Frohhch one when k — » 0. The low-energy part is the delta-functional peak 

V2 \ 



goiuj < 0) 



(At + V2)- 



-S \ to + a- 



V2. 



(5.3) 



(5.4) 



The comparison of o ur numeric results for the low-energy part of the Frohhch polaron (k = 0)spectral density 
for a = 0.05 and Eq. (5.4) demonstrates a perfect agreement (with the accuracy 10^** for the polaron energy and 
Z-factor), whereas our results for the high-energy part (upper panel in Fig. n3) significantly deviate from the analytic 
curve. This is not surprising since for Frohhch polaron the perturbation theory expression is diverging as w — > wq and, 
therefore the perturbation theory breaks down. To test the case when perturbation theory is obviously valid we set 
K = 1 and obtained a perfect agreement for both the low- and high-energy parts of g{uj) (lower panel in Fig. n3). We 
note that the high-energy part of g{uj) is successfully restored by our method despite the fact that the total weight of 
the feature is less than 10^^ for a — 0.05. 

One can note that the main deviation of the actual spectrum of Frohlich polaron from the perturbation theory 
result is the extra broad peak in the actual spectral density at w ^ 3.5. To study this feature we calculated g{i-u) for 
a — 0.5, a = 1, and a = 2 (see Fig. n4h. Note, that the peak is seen for higher values of the interaction constant and 
its weight grows with a. Near the threshold, to ^ I, the spectral density demonstrates the square- root dependence 
~ ^/uj — 1 (see the insert). 

To trace the evolution of the peak at higher values of a we calculated the spectral density for a ~ 4:, a — 6, and 
a = 8 (see Fig. ng). At a = 4 the peak at a; '^ 3.5 already dominates in the spectral density. Moreover, a distinct 
high-energy shoulder appears at a = 4, which transforms into a broad peak at cj ^ 9 in the spectral density for 
a = 6. The spectral density for a = 8 demonstrates further redistribution of the spectral weight between different 
maxima without significant shift of the peak positions. One can also see that there is a high-energy shoulder which 
is, probably, the precursor of another peak which would appear for higher values of the interaction constant. 
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FIG. 14. The spectral density of Frohlich polaron for a = 0.5 (dotted line), a = 1 (solid line) and a = 2 (dashed line), with 
energy counted from the position of the polaron. The initial fragment of the spectral density for a = 1 is shown in the insert. 
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FIG. 15. Evolution of spectral density with a in the cross-over region from intermediate to strong couplings. (The polaron 
ground state peak is shown only for a = 8. Note, that the spectral analysis still resolves it, despite its very small weight 
< 10^"'.) The energy is counted from the position of the polaron. 
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The excited states of the polaron were studied within the frameworks of different approaches |p5|-p9[ by calculating 
optic absorption spectra. The light absorption is associated with the transitions from the polaron ground state (/)o 
with k = and E — Eq to the excited states / with Ef which are characterized by the presence of a finite number 
of real phonons along with the polaron. The optic absorption spectrum at the frequency lu is proportional to the 
transition probability 

PH = 27r^(0o I O I ./)(/ I 6 I cf,o)5{Eo - Ef+u;). (5.5) 

(Here O = Er is the electric dipole interaction, E is the electric field.) 
It was shown in the weak-coupling limit i23,Effl that the optic absorption spectrum has a broad peak with the onset 



separated from the polaron state by the optic phonon energy. Our calculations confirm (see Fig. 14) that there are 
no metastable excited states of the polaron in the weak-coupling regime. 

On the other hand, in the strong-coupling limit the existence of the metastable relaxed excited state (RES), i.e. 
the state where the lattice readapts to the new electronic configuration and the polaron-lattice system is in the local 
minimum of the total energy, was predicted [E^-E9[. This state manifests itself as a sharp peak in the absorption 
spectrum which is located at the frequency equal to the energy difference -Bres — Eq. To check the existence of 
RES one can study the spectral density ( ^.5[ ) since although the matrix elements of transition probability (pM and 
spectral density ( |2.5| ) are different, both functions have to demonstrate sharp peaks at the energies of the metastable 
excited states. From Fig. |l^ we conclude that there is no metastable excited state because the width of the peaks 
is comparable with the excitation energy, i.e. with the distance from the polaron ground state. Moreover, according 
to the strong-coupling approaches |^, the excitation energy of the RES state is proportional to a^, whereas peak 
positions in 17(0;) with respect to Eq do not change with a. 

The variational treatment developed in Ref. pl| suggests that in a certain region of a there may exist two different 
stable states of the polaron [the corresponding equations for variational parameters have two solutions] . Our numeric 
study can shed light on this situation. 

To start with, let us discuss what one could observe would the two states really exist. If at some point a = a* 
there occurs a level crossing so that the ground state switches from one state to another, and the two states differ 
essentially in the number of phonons and/or in the effective mass, one would expect at the point a* a sharp change 
of these quantities. The change should be almost jump-like in the case of small hybridization between the two states, 
and look like a smooth cross-over otherwise. Even in the case of sufficiently strong hybridization, one may distinguish 
between two qualitatively different cases: (i) the case when the level separation is less than ojq [and thus both states 
are stable against decay], and (ii) the case when the upper level is in the continuum, and therefore is unstable. 

Strictly speaking, in the case (ii) one can invoke the second level only in some quantitative sense, since there is no 
qualitative difference between the case (ii) and a situation with only one polaron state. These quantitative features 
could be associated with the non-monotonic behavior of the derivatives (with respect to a) of the effective mass 
and/or the mean number of phonons, and, of course, another peak in the spectral density. 

Our study of the spectral densities shows that the case (i) is not realized because we found that g{0 < uj < 1) — 0. 
Therefore, there is no excited stable state in the energy gap between the ground state energy and incoherent continuum. 
Instead, there are several many-phonon unstable states at energies Ef — Eq ^ 1, ^ 3.5, and ~ 8.5. One can speculate 
that these states reveal themselves in variational approaches and can be mistreated as quasi-stable states of the 
polaron. It should be emphasized, however, that the situation does not resemble that of the level-crossing at all, since 
we do not observe non-monotonic behavior of the derivatives (with respect to a) of the effective mass and/or the 
mean number of phonons. 
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APPENDIX A: UPDATING PROCEDURES 
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1. Updates of class I 

Updating procedures of this class are the simplest. They mimic standard rules of simulating a given distribution 
function Vm ■ In the present case we are dealing with quite a number of variables having different physical meaning: 
external variables {y} include r, N, a, and fc, and internal variables describe the topology of the diagram (index £,m), 
times of electron-phonon vertices and momenta of phonon propagators. From this list of variables follows a set of 
updates simulating multi-dimensional distribution Vm- 

a. Vertex shift in time 

We choose at random any interaction vertex inside the graph (we exclude the diagram closing points which are 
updated separately), and change its time position from t„ to r^ on the interval (ti,T2) between the nearest left- and 
right-neighbor vertices, i.e., ti < r„, t^ < T2. Let the incoming and outgoing electron momenta for the selected vertex 
arc p and p + q. The normalized probability density to find the vertex at time r^ is a simple exponential function 

where AE = E{p) — E{p+q)^ujp depending on whether the updated vertex is the left or right end of the corresponding 
phonon propagator, which allows a trivial solution of the equation 



W{s)ds = r , (A2) 



in the form 



<^n- '-"-"^"'""^» . (A3) 

Here and below r is the random number homogeneously distributed on the unit interval. Since the new variable is 
selected according to the exact probability density the acceptance ratio for this update is unity. 

b. Change of transferred momentum angle 

We choose at random any phonon propagator except those attached to the diagram ends (propagators attached to 
the diagram ends appear in pairs with equal momenta, thus single propagator updates do not apply to them) and 
change its momentum q — > q' so that |q| = |q'|. Let the propagator connects vertices at times ti and T2. Evaluating 
the average electron momentum between these vertices 

J^"p{T)dT 

< P >ri,r,= — , (A4) 

T2 - Tl 

and introducing vector po =< p >ti.t2 +<1j '^6 may write the probability density to find azimuthal and polar angles 
if, 9 between vectors q and po as 

W{ip,9)^sin{0)exp\-^^^^^poqcos9\ , (A5) 



This result is a trivial consequence of the quadratic dispersion law for the bare electron spectrum. Clearly, the new 
azimuthal angle is selected at random {ip — 27rr), and cos 9 is selected according to the simple exponential function 
in complete analogy with Eqs. (|Al|) and ([A^ ) up to trivial change of notations. The acceptance ratio is thus unity. 
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c. Change of transferred momentum modulus 

In this procedure propagators are selected as explained in the previous subsection, but now we change the modulus 
of the transferred momentum while keeping the polar and azimuthal angles between the vectors po and q fixed. The 
probability density now reads 

W^(g)~F(g)g2exp|-^5^[(7-poCos0]2|~exp|-^^^[g-poCOS0]2| , (A6) 

where we have used explicitly the property of the Frohlich model that V{q)q^ is g-independent. By tabulating the 
inverse error-function we ensure fast numerical solution of the equation crrf(z) — r, or z — errf^ (r), and thus 
generation of the new value q = zy^2m/{T2 — ti) + po cos 9 with acceptance ratio unity. 

d. Change of diagram structure 

We select at random any nearest-neighbor pair of vertices inside the graph (again, the diagram closing vertices are 
excluded) and exchange the assignment of the phonon propagators between these vertices. Namely, if the original 
momentum transfer was qi in vertex 1 and q2 in vertex 2, we suggest to change these momenta to q2 and qi 
correspondingly. The acceptance ratio for this procedure depends on whether we are dealing with the left (c = 1) or 
right (c = —1) ends of the phonon propagators 

Jj> _ g--r[-E(p+Ciqi-C2q2)-S(p)-LJo(ci-C2)] fA7) 

where r is the time difference, and p is the electron momentum between the selected vertices. Clearly, this procedure 
effectively changes the topology of bosonic lines while keeping fixed their momenta. 

e. Change of diagram length in time 

This procedure is done in two variants (almost identical to the procedure of shifting the vertex position in time) . 
Consider the case when no artificial potential except the chemical potential is used. In the first variant we select the 
new time difference r between the positions of the right diagram end at its left nearest neighbor vertex according to 
the probability density 

W(r') = ASe-"'^^, (A8) 

where AE — E{p) + NzLUq ~ Mi ^^'^ P is the momentum of the last electron propagator, /i is the chemical potential, 
and Nz is the number of phonon propagators attached to the diagram right end (obviously, Nz is the same for the 
diagram left end). In the second variant we select new time differences between the positions of all nearest- neighbor 
pairs of vertices. For each such a pair the probability density is still given by Eq. ([A^), where in the most general 
case Nz must be understood as the number of phonon propagators which are cut when the diagram is cut anywhere 
between the selected pair of vertices. Notice, that the second variant requires much longer computation time; thus 
if the typical diagram order is very large it must be applied less frequently. In both variants the acceptance ratio is 
unity. 

There is a bottle-neck in the time decay of the electron P-function, which does not allow efficient sampling of both 
long-time and short-time behavior and causes normalization problems at large a, namely, P{t) drops to almost zero 



vaalue at short times, and then climbs back to P ^ 1 before it settles to the asymptotic decay ( 2.15 ). 

There is however a general prescription of how to eliminate such difhculties by using the so-called "guiding function" 
|pO| or fictitious potential renormalization. This method was successfully applied recently |31t| to the problem of 
tunneling transition amplitudes, where one is bound to collect reliable statistics which varies by orders and orders of 
magnitude between different points in time. The idea is to modify the statistics of suggested diagrams by introducing 
the acceptance ratio 

R = Afic(r„cw)Mfic(Toid) , (A9) 

and accordingly multiplying all MC estimators in the time domain by l/AadT), where the fictitious potential j4fic(T) 



is arbitrary. Note, that in Eq. (A9) we are dealing with the external variable r - the diagram length in time. In 
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the present case the best choice would be Afic ~ 1/P(r). We achieve this goal by self-consistently adjusting A^c to 
^/Pmc{t) after a certain large number of updates during the thermalization stage (here Pmc{t) is the statistical 
result for P{t)). After thermalization stage we start collecting new statistics for P{t) and keep Afic fixed. 



/. Change of coupling constant 



Since the diagram weight depends on the coupling constant as a ^ where Np is the number of phonon propagators 
in the diagram, all we need to do is to select new value of a with this power-law probability density. Normalized 
probability density is obtained by restricting allowed values of a to a certain parameter range. Acceptance ratio is 
unity. 



g. Change of external momentum 



Given the average electron momentum of the diagram p =< p >o,r, see Eq. ( A4), with external momentum k, we 
define vector po = k— < p >o,t and write the probability density to select new external momentum k' as 



W^(k')~exp{-^(k'-po)2} 



(AlO) 



As before the new variable is seeded according to this probability density utilizing the tabulated error-function [see 
subsection Ale and Eq. (A6)] and thus is always accepted. 

One note is in order here. Although one is allowed to change the coupling constant and the external momentum 
in a single MC process, it seems more efficient to keep these variables fixed instead of spreading the statistics over 
some range in the (k, a) parameter space. However, the knowledge of the relative weights according to which a 
given diagram contributes to the statistics of various a and k may be utilized in collecting statistics for the finite 
neighborhood of the point (ko,Q;o) used in a given MC simulation. Obviously, reliable results for points other than 
(ko, Qfo) are obtained only provided that for typical diagrams the relative weights are of order unity. As explained in 
the text, this knowledge is also used in deriving estimators for the effective mass and group velocity of the polaron. 



2. Updates of class II 



These updates are in the heart of the method since they change the diagram order. The generic rules 
for constructing them are as follows ||I^. Let the update A transform a diagram Vm{^m,y,xi, . . . ,Xm) into 
T^m+n{^m+my,xi, ■ . ■ ,Xm,Xm+i, ■ ■ ■ ,Xm+n), and. Correspondingly, its counterpart B perform the inverse transfor- 
mation. For n new variables we introduce vector notation: x = {xm+i,Xm+2, ■ ■ ■ ,Xm+n\- The update A involves 
two steps. First, it proposes a change, selecting a new diagram, Vm^n, and a particular value of x, which is seeded 
with a certain normalized distribution function W{x). There are no requirements strictly fixing the form of W{x)^ 
but to render the algorithm most efficient, it is desirable that W{x) be chosen as close as possible to I?m+n(x), i.e., 
to the actual statistical probability density of x in the new diagram. Upon proposing the modification, the update 
is accepted, with probability, Pacc{x), or rejected. The update B, removing variable x, is accepted with probability 
-Prem(a^)). For the pair of complementary updates to be balanced, the following Metropolis- like prescription should be 
fulfilled |l|]: 

" R{x)/W{x), if R{x) < W{x) , 

1 , otherwise , 



-'accl.'^j 



(All) 



^rem\X j 



where 



Wix)/R{x) 
1 , 



,^, PB T), 

R(x) — 

PA 



if R{x) > Wix) 
otherwise , 



i+n (^SrM+ni V^ X\ , . . . , Xm^ X j 
^m \^m 5 y? Xl , . . . , XniJ 



(A12) 



(A13) 



and pj( and pg are the probabilities of selecting updates A and B, which, in principle, may differ. To solve the polaron 
problem and account for any possible diagram it is sufficient to have two pairs of complementary processes of type II 
which are described in detail below. 



23 



a. Adding/removing phonon propagators to the diagram 

Consider the algorithm for the process increasing the number of internal phonon propagators (i.e., excluding those 
attached to the diagram closing points) by one. This update is done in two variants which differ in the probability 
densities according to which the new propagator parameters are suggested. First we select the time position ri for the 
left-hand end of the extra phonon propagator. This is done by choosing at random (with equal probabilities) one of 
the free-electron propagators, and by taking for ri any time (with equal probability density) within this propagator. 
Then we select the transferred momentum and propagator length in time using the distribution function 

H^(q, r) = ^^e-™°(i+9/«'')' (A14) 

Anqo 

where gg/S — wo, i-e., we first seed \q\ according to Wi{q) == l/[qo{^ + q/qo)^] (and isotropic around the point q = 0), 
and then r according to W2((7|r) — ujo{1 + q/qo)^e^'^^°^^^'^^'^°^ . Since the typical length of the phonon propagator in 
time depends on how close is the polaron momentum to the dispersion law end-point, we also use another variant of 
seeding new variables q and t, namely, we factorize the distribution function into W{q, r) = Wi{q' = |q — k|)W3(T) 
(i.e., isotropic around the point k), where W3(r) = fle^'^^^ and ft ^ loq down to 17 ^ O.OIujq close to the end-point. 

We underline, that the above choices are motivated by the physics of the problem, in particular, if the combination 
V^{q)q^ was some power law function of q (e.g., when the interaction vertex is non-singular at small momentum or 
even goes to zero as g -^ 0) one would better have to choose Wi{q ^ 0) oc V^{q)q'^ to ensure that nowhere in the 
accessible parameter region the acceptance ratio (see below) is singular. 

Now the proposing stage is completed, and we are ready to perform accept/reject step, following the above pre- 



scription, Eq. (All). The corresponding function W{x) {x = {Ti,T2,q}) reads (for the first version) 

W{x) = ^^ , (A15) 

where tq is the length of the free-electron propagator, where the point ri is selected. As mentioned already, this form 



of W is by no means the unique one. Apart from the factor pb/pa which will be discussed later, the ratio (A13) is 
now completely defined since 

T^m+n i^m+n , y, Xi, . . . , Xm, x) ^ 2V27ra ^_(T.2_.ri)[c^0+B(<P>xi ,^2 -q)-B(<P>xi,x2 )] (Al6) 

VmiCm^y.Xl,.. . ,Xm) (27r)3 

The algorithm for the process B is to select at random (with equal prob abilit ies) so me p honon propagator, and, if 



it is not attached to the diagram end, with the probabilities given in Eqs. ( A12 ) and ( A15 ) remove it. 

To complete the description of the sub-processes A and B, we should define the ratio pb/pa- It is quite reasonable 
to select creation and annihilation procedures with equal probabilities. At the first glance it might seem that this 
immediately leads to pb/pa — Ij but this is not true. The point is that when we select an electron propagator for 
placing the point ti , we have Ng equal chances, where N^. is the number of free-electron propagators in the diagram 



being modified [denominator of Eq. (A_13) ], and when we select a phonon propagator for removing, we have Nph equal 
chances, where Npf^ is the number of phonon propagators in the diagram from which we try to remove the propagator 



[numerator of Eq. ( A13| ) ]. These N,, and Np^ are straightforwardly related to each other: 

Nph = (N, + l)/2 . (A17) 



(A18) 



We thus get 

PB _ 2JVe ^ ^Nph - 1 

PA Ne + l Nph 

(Note a misprint in Ref. 0, where the r.h.s. of Eq. (10) gives pa/pb instead oi pb/pa)- 

b. Adding/removing a pair of phonon propagators attached to diagram ends 

This update is done in close analogy with the previous one, except minor changes to which we proceed now. 
First, we select time positions ri and T2 for the left- and right-end propagators according to the probability densities 
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Wi{ti) — fie ^^^ and Wr{T2) = fie ^^"^ ^^\ In the first variant of the update fl = luq and in the second fl <C wq- 
Also, the phonon momentum is suggested using the same distribution Wi{q' = q) or Wi{q' = |q— k|). We thus have 

^2 -n{T + Ti-T2) 

Wix) = — — r , (A19) 

and 



Vyn+n{£,ni+n,y, Xi,...,Xm,x) _ 2^/27ra -^„(r+Ti-T2) 
Vjn{£.m,y,Xi,...,Xra) (27r)3 

g-[£;(<p>o,xi-q)-B(<p>o,ri)]ri-[B(<p>^2,--q)--E(<P>x2,x)](r-r2) ^^ ^ ^^ 



g-[£;(<p>0,x2-q)--E(<P>0,^2)l-^2-[B(<P>x2.xi-2q)-B(<p>,2,^i)](Ti-r2)-[£(<p>,,.,-q)-£;(<p>,i,,)](r-ri) ^ > T2 ' (^^'^-' 

The algorithm for the inverse procedure is to select at random (with equal probabilities ) a p air of propagators from 
the list of pairs attached to the diagram end, and with the probabilities given in Eqs. ( |A12| ) and ( A19 ) remove it. 



Since we select procedures inserting and removing pairs of propagators with equal probabilities, we have 

PB _ 1 
PA~ N^ + l' 



(A21) 



APPENDIX B: METHOD OF SPECTRAL ANALYSIS 

1. General background and outline of the method 

The problem of restoring positive definite spectral density function p{(jj) from known imaginary-time Green function 
G{t) is the problem of solving linear first- type Fredholm equation 



e-'''^ p{uj)duj = G{t) , (Bl) 

where the domain of definition of the functions G{t) and p{w) is [0,oo]. The normalization of the Green function 
G(0) = 1 implies the additional constraint 

/•CXD 

p{uj)duj = 1 . (B2) 







In a realistic situation the Green function is known at a discrete set of times {Ti},i = 1, ...,N with some statistic 



errors at each time point. As is well known, in this case the problem of solving Eq. (Bl) belongs to the class of 



ill-posed problems. The characteristic feature of the ill-posed problem is that the solution of equation (Bl 



IS not 



unique even when statistic errors are absent, as there is an infinite number of unknown functions p{uj) satisfying (Bl) 



In the case of finite statistic errors one may face a situation when the solution of Eq. (Bl) under the constraint (|B2| ) 
does not exist at all. Therefore, it is natural to formulate the problem as to find an approximate solution Pmin{^) 
which reproduces G at a finite set of times with smallest deviation -Dmin- The definition of the measure of deviation 
depends on the method used, and the value of minimal deviation Dmin is determined by the magnitude of statistic 
errors. 

There are two fundamental difficulties that are inherent to the spectral analysis. The first one is the well-known saw- 
tooth instability of the linear Fredholm equation of the first type - an approximate solution p{uj) does not reproduce 
the true solution p(w) even if p{uj) generates the Green function 

/•oo 

G{t) = / e-'^^ p{uj)duj (B3) 

Jq 

which reproduces G{t) with any preassigned accuracy. This difficulty is treated usually by the rcgularization method 
that smoothes the saw-tooth noise of approximate solution p{lu). The idea of the rcgularization method is to introduce 



some nonlinearity into Eq. (Bl) that imposes constraints on the derivatives of p{uj). There are two main drawbacks of 



this method. First, regularization method is unable to restore the spectral density which has sharp features. Second, 
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due to a distortion of the initial equation by additional regularization terms the approximate solution reproduces the 
function G{t) with relatively high deviation D 3> i'min- Hence, the information from the most representative region 
of the deviations D ~ Dmin is lost. 



The second difficulty inherent to the problem of solving equation (Bl) is that any representation of p(w) by a 
preassigned discrete set {p{ujf)}, f — I, ...,M is the source of uncontrollable systematic errors. For one thing, if the 
function p{uj) contains a sharp feature with a significant weight at some w', which does not match the discrete set 
{ujf}, this feature cannot be reproduced properly and, therefore, the rest of the spectral density can be distorted 
beyond recognition. Note, that all iteration methods as well as the methods based on solving the nonlinear system of 
equations use preassigned discretisation of the w-space. 



We present a method of solving equation (Bl) that avoids distortion of equation by nonlinear terms and, thus, probes 
the most representative interval of deviations. Besides, the method does not suffer from systematic errors as it does 
not involve preassigned discretisation of the w-space. The idea of the method is to generate by a stochastic procedure 
a (large enough) set of M positive definite statistically independent approximate solutions {pj{uj)},j = 1, ..., M with 



deviation measures Dj ^^ Dmin- And then, taking advantage of the linearity of Eq. (Bl), choose the final solution as 
the average 



M 



p{uj) ^ M-^^p,{u:) . (B4) 

The reason is that while the particular solution Pj{lo) possesses the saw-tooth instability, the stochastic character 
of the procedure of particular solution generation should lead to averaging out the saw-tooth noise. Note, that the 



condition Pj{uj) > and constraint ( p32| ) substantially enhance the convergence of the averaging (B4) 
The method of generation of a particular solution is based on the optimization of the deviation 



D[p] = 



Tniax 



G{t)~G{t) G-\T)dT . (B5) 



is the maximal r up to which G{t) is known. The weight function G~^{t) is to efficiently utilize information 
from the whole range [0,rinax], even in the case when the function G{t) decreases by orders of magnitude with r. 
Note, that we use weight function G~^{t) rather than G^^{t) to avoid feedback instabilities in generation p(w). 

Our optimization procedure does not involve preassigned fragmentation of the cj-space. The number of parameters 
used for parametrization of the spectral density function p{lo) is being varied during optimization process, so that 
any spectral function can, in principle, be reproduced within any preassigned accuracy. The process of generating a 
particular solution involves a random choice of the initial-configuration parameters and subsequent optimization of 
the deviation by changing the parameter values, as well as the number of the parameters. The maximal number of 
continuous parameters and the number of particular solutions M are limited only by the computer performance. 

2. Configuration and method of getting independent solution 

We parametrize p as a sum 

K 

t=i 

of rectangulars {Pt} — {ht^wt, Ct} 

^iP^}(^^-[0 , otherwise. (^^) 

determined by height hf > 0, width wt > 0, and center Cf > 0. 
A configuration 

C = {{Pt},t = l,...,K} (B8) 

with the constraint 
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K 

J2 ^i^t = 1 (B9) 



defines, according to Eqs. (B6, B3), the fmiction G{t) at any time point 



1 , r-0, 

Gc{r) = { 2^_i g ^^^_,,, gi^i^(^^^/2) , r ^ . ^^^°^ 

To express the deviation ( |B5| ) as an analytic function of the values of G and G at the set of times {Ti\, i — \, ..., N 
[where the function G{t) is known], we use linea r in terpolation between closest points. 



Note, that the specific type of the fu nctio ns (B7) is not crucial for the general features of the method although 
simple form of analytic expressions (B9, |B10|) is of value for fast performance. 

The procedure of obtaining a particular solution Pj{lo) consists of randomly generating some initial configuration 
C'"'* followed by nondeterministic sequence of configuration changes until deviation satisfies the condition 

D[Cf^] <Du-^ i^min (Bll) 

{Du is the upper limiting deviation.) for final configuration C^". The nondeterministic character of configuration 
changes is achieved by random selection of various elementary updates. 

3. General features of elementary updates 

By elementary update we mean a random change of the configuration, which is either accepted or rejected in 
accordance with certain rules. There are two classes of elementary updates. The updates of the class I do not 
alter the number of rectangulars, K, changing only the values of the parameters from a randomly chosen set {Pt\- 
The updates of the class II either add a new rectangular with randomly chosen parameters {/ii<-+i,i(Jif+i, ca'+i}, or 
remove stochastically chosen rectangular t from the configuration. If a proposed change violates constraint (B10|) 



(e.g., a change of Wt or ht, or any update of the class II), then the necessary change of some other parameter set {Pf} 
is simultaneously proposed, to satisfy the requirement of the constraint. 

The updats should keep parameters of a new configuration within domain of definition of configuration C. Formally, 



the domains of definition of configuration (B8) are S/j^ = [0,oo], 'B.c^ — [0,cxd], 2^,^ — [0, 2ct], and S^ G [l,oo]. 
However, for the sake of faster convergence, we reduce domains of definition. 

As there is no general a priori prescription for choosing reduced domains of definition, the rule of thumb is to 
start with maximal domains and then, after some rough solution is found, reduce the domains to reasonable values 
suggested by this solution. In particular, since the probability to propose a change of any parameter of configuration 
is proportional to K~^ , it is natural to restrict maximal number of rectangulars 'E.k € [1, -?^max] by some large number 
^max- To forbid rectangulars with extremely small weight, which contribution to G{t) is less than statistic errors of 
G{t), one can impose the constraint htWt € [S'min, 1], with S'min ^ -K'max- When there is some preliminary knowledge 
that overwhelming majority of integral weight of the spectral function p{ijj) is in a range [wmin, Wmax], one can restrict 
the domain of definition of the parameter ct by Scj = [wmin, Wmax]- Then, to reduce the phase space one can choose 

'^ht = [Vin,C)0] and 2^,^ = [■U;inin,min{2(ct - CJmin), 2(tJniax - Ct)}]. 

While the initial configuration, the update type, and the parameter to be altered are chosen stochastically, the 
variation of the values of the parameters relevant to the update is optimized to maximize the decrease of T). Each 
elementary update of our optimization procedure (even that of the class II) is organized as a proposal to change some 
continuous parameter ^ by randomly generated 5S, in a way that the new value belongs to S^. Although proposals 
with smaller values of 5^ are accepted with higher probability it is important, for the sake of better convergence, to 
propose sometimes changes 5^ that probe the whole domain of definition Sj. To probe all scales of (5^ G [(^^min, <^6nax] 
we generate 5£^ with the probability density function V ~ (max(|(5^i„in|, |<JCmax|)/|'^CI)^> where 7 :^ 1. 

Calculating the deviation measures D{£_), D{^ + 6^), D{^ + 6^/2), and searching for the minimum of the parabolic 
interpolation, we find an optimal value of the parameter change 

S^opt = -B/2A, (B12) 

where 
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A = 2{D{i + 50 - 2D{^ + 5i/2) + D{0)m- 



(B13) 



and 



B = (4i?(e + 6i/2) - D{i + 5i) - 3D{0)6^. 



(B14) 



In the case A> and ^opt G ^^ we adopt as the update proposal S£_ one of the values S£_, S^/2, or S£_opt for which the 
deviation measure D{^ + S^ is the smallest. Otherwise, if the parabola minimum is outside S^, one has to compare 
only deviations for 6^ and S(_/2. 



4. Global updates 



The updating strategy has to provide efficient minimization of the deviation measure until criterion (Bll) is satisfied. 
It is highly inefficient to accept only those proposals that lead to the decrease of deviation, since, in a general case, 
there is an enormous number of deviation local minima Z?ioc[C] > -D«. As we observed it in practice, these multiple 
minima drastically slow down (or even freeze) the process. 

To optimize escape from a local minimum, one has to provide a possibility of reaching a new local minimum with 
lower deviation through a sequence of less optimal configurations. It might seem that the most natural way of doing 
this would be to accept sometimes (with low enough probability) the updates leading to the increase of the deviation. 
However, this simple strategy turns out to be impractical. The reason is that the density of configurations per interval 
of deviation sharply increases with D. So that the acceptance probability for a deviation-increasing update should be 
fine-tuned to the value of D. Otherwise, the optimization process will be either non-convergent, or ineffective [if the 
acceptance probability is, correspondingly, either too large, or too small in some region of D]. 

A way out of the situation is to perform some sequence of T temporary elementary updates of a configuration C(0) 

C(0) ^ C(l) ^ ... ^ Cir) -^ C{r + 1) ^ ... ^ C(r) , (B15) 

where the proposal to update the configuration C{r) -^ C{r + 1) is (temporary) accepted with the probability 



Vr 



-tr+1 



f{D[C{r)]/D[C{r + l)]) 



D[C{r + l)]<D[Cir)], 
D[C{r + l)] >D[C{r)] . 



(B16) 



( Func tion / satisfies boundary conditions /(O) = and /(I) — 1.) Then we choose out of the configurations {C{r)} 
( B15D the one with minimal deviation and, if it is different from C(0), declare it to be the result of the global update, 
or, if this configuration turns out to be just C(0), reject the update. 
We choose the function / in the form 



/(^) 



^l+d 



(d > 0) , 



(B17) 



which leads to comparatively high probabilities to accept small increases of deviation measures and hampers significant 
enlargements of deviation. Empirically, we found out that the global update procedure is most effective if one keeps 
parameter d = di ^ at the first Ti steps of sequence (B15) (to leave local minimum) and then changes this parameter 
to a value d = ^2 3> 1 for the last T — Ti elementary updates (to decrease the deviation measure) . In our algorithm 
the values T £ [l,T'max], Ti e [1,T], di G [0, 1], and d2 G [l,dmax] were stochastically chosen for each global update 



5. Final solution and refinement 



After a set of A/ configurations 



{Cf,j = l,...,M} 



(B18) 



that satisfy the criterion (Bll) is produced, the solution (B4) is obtained by summing up the rectangulars ( B7|jB18 ). 

We, however, employ a more elaborated procedure, which we call refinement. Namely, we use the set ( p318D as a 

source of M^cf new independent starting configurations for further optimization. These starting configurations are 



generated as a linear combinations of randomly chosen members of the set (B18) with stochastic weight coefficients. 
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Then, the refined final solution is represented as the average (B4) of Afrcf particular solutions resulting from the 
optimization procedure. 

The main advantage of such a trick is that initial configurations for optimization procedure now satisfy the criterion 



(Bll) from the very beginning and, thus, upper limiting deviation D^ can be considerably reduced. Moreover, as 
any linear combination of sufficiently large number R of randomly chosen parent configurations IC^",?? — 1, ...,R\ 
smoothes the saw-tooth noise, the deviation of a summary configuration C^^ is normally lower than that of each 
additive one. 

6. Elementary updates of class I 

(A) Shift of rectangular. Change the center ct of a randomly chosen rectangular t. The continuous parameter for 



optimization (B12-B14) is ^ = c* which is restricted by domain of definition Scj — [wmin + Wt/'^, ^max ~ Wt/2]. 

(B) Change of width without change of weight. Alter the width wt of a randomly chosen rectangular t without 
change of the rectangular weight hfWt — const and center ct. The continuous parameter for optimization is ^ ~ Wt 
which is restricted by St„j = [wmin,niin{2(ci - w,„i„), 2(wmax - q)}]- 

(C) Change of weight of two rectangulars. Change the heights of two rectangulars t and t' (where i is a randomly 
chosen and t' is either randomly chosen or closest to t rectangular) without change of widths of both rectangulars. 
Continuous parameter for optimization is the variation of the rectangular t height S, — ht. To restrict the weights of 



chosen rectangulars to [5'min, 1] and preserve the total normalization ( [B2|) this update suggests to change ht -^ ht + S£, 
and hf -^ ht' — d^Wt' /wt with d£, confined to the interval 

Smin/Wf - ht < S^ < {hf - Si-nin/wt')wt/wt' (B19) 



7. Elementary updates of class II 

(D) Adding a new rectangular. To add a new rectangular one has to generate some new set {Pncw} — 
{hnew,Wnow,Cncw} and reduce the weight of some other rectangular t (either randomly chosen or closest) in order 



to keep the normalization condition (B2). The reduction of the rectangular weight t is obtained by decreasing its 
height ht. 

The center of the new rectangular is selected at random according to 

Cnow = (^min + Wniin/2) + (o^max - t^min " Wmin)'' (B20) 

As soon as the value Cnew is generated, the maximal possible width of a new rectangular is given by 

wJJJ,'^ = 2 min(cJniax - Cnow, C„cw " Wmin)- (B21) 

Continuous parameter for optimization 6^ = h^cwWucw is generated to keep weights of both new rectangular and 
rectangular t larger than 5min 

S^ - 5min + r{htWt ~ 5min) (B22) 

Then, the value of the new rectangular height ft-ncw for given S^ is generated to keep the width of new rectangular 
within the limits [waiin,w^c^] 

hnc^ = '5e/<ew + riS^Wmin " '5C/<ew ) ■ (B23) 

(E) Removing a rectangular. To remove some randomly chosen rectangular t, we enlarge the height hf of some 



another (either randomly chosen or closest) rectangular t' according to condition (B2). Since such procedure does 
not involve continuous parameter for optimization, we unite removing of rectangular t with the shift procedure (A) 
of the rectangular t' . Then, the proposal is the configuration with the smallest deviation measure. 

(F) Splitting a rectangular. This update cuts some rectangular t into two rectangulars with the same heights ht 
and widths Wncwi = Wmin + i'{u!t — Waiin) and Wnow2 = Wt — lUncwi- Sincc removing a rectangular t and adding of 
two new glued rectangulars does not change the spectral function we introduce the continuous parameter 5£^ which 
describes the shift of the center of a new rectangular with the smallest weight. Second rectangular is shifted into 
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FIG. 16. Model spectral density (dashed line) and the result of spectral analysis (solid line). The position of the delta- function 
is shown only in the lower panel. 



opposite direction to keep the center of gravity of two rectangulars unaltered. The domain of definition S^ obviously 
follows from the parameters of the new rectangulars. 

(G) Gluing rectangulars. This update glue two (either randomly chosen or closest) rectangulars t and t' into single 



new rectangular with the weight h„ 



vW„ 



hfWt + Wt'hf and width Wn 



{wt + Wt')/2. The initial center of 



the new rectangular Cncw corresponds to the center of gravity of rectangulars t and t'. We introduce a continuous 
parameter by simultaneously shifting the new rectangular. 



8. Tests 



To check the accuracy of our approach, we tested it for the spectral density distribution that spreads over large range 
of frequencies and simultaneously possesses fine structure in low-frequency region. The test spectrum was modeled 
as the sum of the delta-function with the energy es = 0.03 and the weight Zs = 0.07, and continuous high-frequency 
spectral density which starts at the threshold £th = 0.04. The continuous part of the spectrum pcon was modeled by 
the function [In fact, this functional form is predicted by the Pitaevskii theory for p{uj) near the end point.] 



.H = 



Zs^/lo - eth 



2vr^egap[(w - £th) + Cgap] 



(B24) 



(here Cgap = £th — e^ is a microgap) in the range ui G [sth, 0.566] and by a triangule at higher frequencies (see the 
dashed line in the upper panel of Fig. M). The Green function G{t) was calculated from the model spectral density 
in the n^ax = 300 points r^ = Tmax* /'^max in the time range from zero to r^ax = 1000. The restored spectral 
density reproduces both gross features of high-frequency part (upper panel in Fig. [iq ) and the fine structure at small 
frequencies (lower panel of Fig. Hq). The energy and the weight of the delta- function was restored with the accuracy 
lO"**. The final solution was obtained by the averaging (B4) of M = 1100 particular solutions. 
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FIG. 17. The model spectrum (dashed hnes) and results of spectral analysis (solid lines) for rj — 10 ^ (upper panel) and 



T] — 10 (lower panel). 
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To evaluate the precision which characterizes how typical particular solution Pj{uj) (see (|B3[)) reproduces the Green 
function G{t), we introduce the maximal relative deviation 



77 = max 



G{n) - Gin) 
Gin) 



ie [l,?^max] , (B25) 



which typical value is 77 = 10^* for a particular solution of the spectrum in Fig. y^. 

Since the Green function which is obtained from Monte Carlo calculations contains some statistic errors at each 
time point, the minimal value of parameter rj is limited by the quality of the calculated Green function G{t). To study 
the influence of the (uncorrelated) statistic errors we studied the stability of the method against stochastic noise 

G{n) -^ G(t,)(1 + Tin) , « = 1, J^max , (B26) 

introduced by random numbers r^ G [0,1]. It is seen that the method restores the gross features of the spectrum 
(position and width) even for rather roughly calculated Green function, with rj = 10^^ (upper panel in Fig. |l7[ ), 
whereas the precision 77 — 10^^ is sufficient to resolve the lineshape (lower panel in Fig. O)). 
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